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Reconstruction of the current flow in a vessel system 



The invention relates to a data processing device and a method of 
reconstructing the current flow in a vessel system, and to an assembly that comprises a data 
processing device of this kind. 



5 

Known from US 2002/0041654 A 1 is an X-ray device with which the 
propagation of an injection of contrast medium in the vessel system of a patient can be 
observed. The user hereby has to move the X-ray apparatus during the propagation of the 
contrast medium in order to maintain the propagation front in the image. In order to be able 
10 to execute this movement process without a radiation loading on the patient, a simulated 
contrast-medium propagation is shown on a monitor using stored pictures. Modeling of the 
actual current flow is not provided with this system. 

Against this background, it was an object of the present invention to provide 
means for reconstructing the current flow in a vessel system, which means enable modeling 
15 of the actual circumstances as precisely as possible. 

This object is achieved by a data processing device with the features as 
claimed in claim 1, by an assembly with the features as claimed in claim 8 and by a method 
with the features as claimed in claim 10. Advantageous embodiments are contained in the 
dependent claims. 

20 The data processing device in accordance with the invention serves for the 

reconstruction of the current flow in a vessel system. The vessel system may be, in particular, 
a section of the blood vessels that is of interest, e.g. the system of coronary vessels, without, 
however, the invention being restricted to this application example. The data processing 
device comprises a memory containing measurement data describing an observed progressive 

25 propagation of a (first) medium in the said vessel system. This medium may, in principle, be 
any object that executes a progressive movement through the vessel system. The medium 
may be, in particular, a fluid flowing within the vessel system, such as an injection of a 
contrast medium. The propagation of the medium may be described by means of the 
measured progression over time of the position of the medium, for example of the position of 
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a contrast-medium front or a contrast-medium concentration maximum. This data processing 
device is equipped to reconstruct, from the said measurement data, a modeled propagation, or 
"model propagation'*, of a second medium within the vessel system. The second medium 
may be identical with the first, or may differ from it (e.g. contrast medium vs. blood). In 
5 reconstructing the "model propagation", the following two conditions are to be fulfilled for 
the vessel system under observation: 

- firstly, the difference between the observed propagation and the model 
propagation is to be minimal. In other words, the model propagation sought should follow the 
measurement data as closely as possible. 

10 - secondly, the model propagation is to be monotonously progressive in 

relation to the vessel system. This means that, in accordance with the model propagation, a 
medium can move through the vessel system only "forwards" in one direction. 

With the data processing device described, it is possible to reconstruct a model 
consistently and in a stable manner from measurement data concerning the propagation of a 

15 (first) medium in a vessel system. This is accomplished by matching a model function to the 
measurement data globally, i.e. taking into account the entire vessel system under 
observation, wherein the necessary boundary condition of monotonous propagation is 
simultaneously observed. If, for example, the propagation of a contrast-medium injection 
within the vessel system of a patient is measured with X-ray images, strongly divergent data 

20 arises in practice as a result of numerous interference effects. It may be, in particular, hereby 
that no monotonous, i.e. constantly progressing, propagation of the injection is described by 
the data, but, rather, that jumps and/or backflow are occurring. By contrast, monotony is 
strictly observed when calculating the model propagation for the entire vessel system, as a 
result of which, on the one hand, a compulsory boundary condition is fulfilled and, on the 

25 other, a significant stabilization of the current reconstruction is achieved. 

In accordance with a first further embodiment of the data processing device, 
the reconstruction is executed in such a way that the model propagation additionally has as 
smooth as possible a progression. If the model propagation is described generally as a 
function f(x) of a variable x, the "smooth progression" may be defined, in particular, to the 

30 effect that the sum of the second derivative, |f *|, is to be minimal on average. If f(x) is, for 
example, the position of a propagation front at time x, then the smoothness clearly 
corresponds to a propagation that is as low-acceleration as possible. 

In accordance with a further preferred embodiment of the data processing 
device, the memory contains, as measurement data, bolus arrival times mj. The vessel system 
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is hereby logically divided into individual sections, which are identified by the indices 
i=l,...N, and the bolus arrival time mi is defined as the time determined in a measurement 
which a (first) medium requires, starting from a predetermined starting point, to reach the 
vessel section with index i. On the one hand, the discretization of the problem on which this 
5 embodiment is based is well suited to handling in a digital data processing device, on the 
other, it also accords with the fact that the observed propagation of a medium is generally 
sampled by means of measurements at discrete moments. 

As a further development of the embodiment described above, the data 
processing device is further equipped to calculate model bolus arrival times U for the vessel 
10 sections i in such a way that the following two conditions are fulfilled: 

A, Vi = WV-l (1) 

and 

E = -t\ is minimal. (2a) 

Equation (1) expresses the requirement for a monotonously progressing 
1 5 propagation, whilst (2a) requires, globally for the vessel system, that the model bolus arrival 
times ti should follow, as closely as possible, the observed bolus arrival times mj. The values 
p(i) each hereby reflect the index of the vessel section which, viewed in the direction of flow, 
is located in front of the vessel section with index i. 

In addition, the above-described data processing device may also be 
20 supplemented to the effect that the following requirement is additionally taken into account: 

_ " 

E m =Z-s f i is minimal (2b) 

wherein the quantity I contains the indices of all "internal" vessel sections 

ft 

with a predecessor and a successor, and /, is the discrete approximation of the second 
derivative in vessel section i. The condition (2b) expresses the requirement for as smooth as 
25 possible a progression of the model bolus arrival times tj. With the combination of the 
requirements (2a) and (2b), a weighting factor will generally have to be assigned, the 
defining of which enables the relative significance of the two criteria to be weighted as 
desired. 

In the case of a further development of the embodiment described above, the 
30 data processing device is equipped to calculate the model bolus arrival times t\ using the 
algorithm of linear programming. In the formulation according to equations (1), (2a) and, 
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where applicable, (2b), the model bolus arrival time sought may be regarded as the solution 
to a linear optimization problem, so recourse may be had to algorithms of linear 
programming developed for the purpose, such as the Simplex method. 

In particular, a three-dimensional model of the vessel system in which the 

5 current flow is to be reconstructed may be stored in the memory. Three-dimensional vessel 
models are increasingly encountered within the context of medical applications, and are of 
great value in the interpretation of data by the administering doctor. 

In accordance with another further embodiment of the data processing device, 
it may be equipped to take account of physiological boundary conditions in the 

10 reconstruction of the model propagation. For example, during an injection of contrast 

medium, account may be taken of the preservation of the injected contrast-medium mass on 
its propagation in the vessel system if the measurement data used to describe the propagation, 
or the model propagation sought, makes provision for the mass value to be taken into 
account. By taking account of physiological boundary conditions that necessarily have to be 

15 fulfilled, the calculated model propagation may be further stabilized against interference 
effects. 

The data processing device is further preferably coupled with a display device, 
such as a monitor, in order that the calculated model propagation may be graphically 
represented there. In this manner, the results of the calculation may be made available 

20 visually to, for example, an administering doctor in a catheter investigation. 

The invention further relates to an assembly for observation of the current 
flow in a vessel system, wherein the assembly comprises an image-generating device for 
generating images of the vessel system, and wherein measurement data describing the 
progressive propagation of a medium in the vessel system may be obtained from the said 

25 images. The assembly further comprises a data processing device of the type explained above 
for reconstructing the current flow in the vessel system. This means that the data processing 
device comprises a memory with measurement data, which describes the observed 
progressive propagation of a (first) medium, and which has been obtained by the image- 
generating device. The data processing device is further equipped to calculate, globally for 

30 the vessel system under observation, a model propagation of the (second) medium, which 
minimizes the difference between the observed and the modeled propagation, and 
simultaneously describes a monotonous propagation. 

The image-generating device may be, in particular, an X-ray apparatus, 
preferably a rotating X-ray device for generating three-dimensional X-ray images of a body. 
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Equipment of this kind is generally available as a matter of course in catheter laboratories 
and can therefore be used without problems for a determination of the propagation of an 
injection of contrast medium and thereby of the flow characteristics in the vessel system of a 
patient. 

5 The invention further relates to a method of reconstructing the current flow in 

a vessel system, comprising the following steps: 

a) Obtaining measurement data describing an observed progressive 
propagation of a medium in the vessel system. 

b) Reconstructing a model propagation of a medium in the vessel system in 
10 such a way that: 

- the difference between the observed propagation and the model propagation 
is minimal, and 

- the model propagation is monotonously progressive. 

The method implements, in a general form, the steps that can be executed with 
15 the above-described data processing device or assembly. To explain the method and the 
variants that are possible, you are therefore referred to the above description. 



The invention will be further described with reference to examples of 
20 embodiments shown in the drawings, to which, however, the invention is not restricted. 

Fig. 1 shows, schematically, the components of an X-ray assembly in 
accordance with the invention for reconstructing the current flow in a vessel system. 

Fig. 2 shows the definition of variables and indices in respect of the three 
dimensional model of a vessel tree. 
25 Fig. 3 shows some of the equations required for the linear programming. 



Fig. 1 shows, as a typical example of the use of the present invention, a 
medical X-ray device. This comprises an X-ray apparatus 1, which is affixed to a C-arm and 
30 can rotate about a patient 2 in order to produce X-ray images. The recorded data is passed on 
to a data processing device 3 for further processing. 

The X-ray apparatus 1 may, in particular, be equipped to generate three- 
dimensional reconstructions of a body volume in order that, for example, the spatial structure 
of a vessel system may be investigated. Further, the attempt has been made to reconstruct 
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sequences over time, such as, in particular, the flow or propagation of a fluid in the vessel 
system (see H. Schmitt, M. Grass, V. Rasche, O. Schramm, S. Haehnel, K. Sator: "An X-ray- 
based method for the determination of the contrast agent propagation in 3-D vessel 
structures", IEEE Trans. Med. Imag., 21(3):251-262, 2002). These kinds of "4D images" 
5 encounter considerable technical problems, for example the necessity for high image rates to 
record the current flow, the stability of the volume segmentation (i.e. assignment of voxels to 
objects) and movement correction. In addition, there is the fundamental problem of a good 
flow reconstruction for an entire vessel system. In this regard, according to the prior art, 
bolus arrival times (BAT) are determined from the X-ray images for the vessel tree, and 

10 assigned to the vessel tree. However, the data thereby obtained shows a severe scatter, which 
is due inter alia to noise in the X-ray images used, inaccuracies during segmentation, 
incomplete movement correction or imprecise determination of the bolus arrival times. 

Against this background, a method for a stable, reliable reconstruction of the 
flow characteristics in a vessel system will be explained below with reference to Figs. 1 to 3. 

15 The method begins with a three-dimensional reconstruction of the vessel tree 6, based on 
rotatory 3D X-ray images. Details of this reconstruction may be found in the relevant 
literature (e.g. M. Grass, R. Koppe, E. Klotz, R. Proksa, M.H. Kuhn, H. Aerts, J. Op de Beek, 
R. Kemkers: "3D reconstruction of high contrast objects using C-arm image intensifier 
projection data", Comput. Med. Imaging Graph., 23:31 1-321, 1999; M. Grass, R. Koppe, E, 

20 Klotz, J. Op de Beek, R. Kemkers: "3D reconstruction and imaging based on C-arm 

systems", Med. Biol. Eng. Comp., 37(2): 1520- 1521, 1999). Following segmentation of the 
vessels, the structure of the vessel tree 6 is known. 

In accordance with Fig. 2, the vessel tree 6 may be divided into individual 
sections or clusters of volume pixels (voxels), which are identified by an index i=l,...N. Per 

25 definitionem, the greatest index N is assigned to the starting point of the vessel tree. Further, 
the function p(i) in each case reflects the index of the vessel section located directly in front 
of the vessel section with index i, viewed in the direction of flow. This function is necessary, 
in particular in the case of junction points in the vessel tree, because a consecutive numbering 
of both vessel branches in steps of units is not possible here. 

30 Once the vessel tree has been established using a first series of angiographic 

X-ray images, a second image series is produced, in which the X-ray apparatus 1 is 
preferably stationary, i.e. does not rotate around the patient 2. 

During the second image series, a contrast medium is injected - e.g. via a 
catheter - into the vessel system of the patient 2, and the resultant propagation of the contrast 
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medium in the vessel system is followed on the generated X-ray projections. Using image 
processing methods, for every section i of the vessel tree, the associated bolus arrival time m* 
may be determined from the projection images. This is the time required by the contrast 
medium from the moment of injection until the vessel section i is reached. The bolus arrival 

5 times may serve for describing the flow characteristics in the vessels, and from them may be 
calculated, for example, the local flow rates. As already mentioned, however, the measured 
bolus arrival times m\ are subject to a high degree of noise as a result of numerous 
interference effects. They therefore frequently infringe fundamental physiological boundary 
conditions, such as the monotony of the progressive flow propagation. This means that a 

10 vessel section located further back in the direction of flow may have a shorter measured 
bolus arrival time mi than a section location in front of it. 

To eliminate this problem and to generate a stable reconstruction of the flow 
characteristics in the vessel system, the calculation of a model propagation is proposed in the 
form of model bolus arrival times ti which, in the context of a global optimization, fulfill the 

15 monotony condition and are simultaneously as close as possible to the measured data. With 
the nomenclature introduced above, the monotony condition for the model bolus arrival times 
ti sought may be expressed by the following equation: 

A,=',-' p(O >0 V,=1,..JV-1 (1) 

Observing this boundary condition, the cost function: 

20 * = XK-',| (2a) 

which quantifies the discrepancy between the measured and modeled bolus 
arrival times, as calculated over the entire vessel progression, is now to be minimized. The 
calculation of the sought variables ti may hereby take place using linear programming (see 
D.G. Lunenberger: "Linear and Nonlinear Programming", Addison-Wesley, Reading, 2 nd 
25 edition, 1984; W.H. Press, S.A. Tuekolsky, W.T. Vetterling, B.P. Flannery: "Numerical 
Recipes in C", chapter 10.8, University Press, Cambridge, 2 nd edition, 1992). With the 
normal form of this technology, the aim is to minimize the function: 



F = ca (3) 
while observing the boundary conditions: 

30 55=6 (4) 

a,>0 (5) 
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wherein b , c and S are given vectors, or a given matrix, and the vector a 
represents the solution sought. 

To be able to use the technology of linear programming for the present 
problem, two N-dimensional auxiliary vectors s* and e" and an N-dimensional auxiliary 

vector A are to be introduced, and the following vectors a , b , c as well as the matrix S are 
to be defined: 



a = 



r 



£ 

A 



f in \ 



'6 



C — 



N 



10 



V P ®{N-\)xN^{N-\)xN ~ l(^_|)x(A^-l) J 



(6) 



(7) 



(8) 



(9) 



15 



Here, l w and 0 M designate M-dimensional vectors, which contain a one or a 
zero in each component; Imxm is a unit matrix of the variable MxM; Omxl is a matrix with M 
rows and L columns, whose components are all zero; and P is the (N-l)xN coupling matrix, 
which is defined by: 

1 ifi = j 

P u =|-i»/ P(0 = J 00) 
0 otherwise 

Since p(i) * i , the definition of P is uncontested. 
Using the variables defined above, the function: 



ft , . 



(11) 



20 



is to be minimized by linear programming, observing the boundary conditions: 
t f + e* - s7 = m i for i = 1... JV (12) 
f,-/„ (f) -A,= 0 /or i = l,.JV-l (13) 
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/,>0 (14) 
<>0 (15) 
e;>0 (16) 
A,>0 (17) 
5 Equation (11), together with equations (12), (14)-(1 7), implies that, for each 

index i=l,. . .N, either e* or ej has to be zero. If, for example, eT > e* > 0 , with the formula 

ef := sj f| + =0a smaller value would be obtained in equation (11). Since in addition 

=m,-/, (18) 

the minimization of the variable F in equation (1 1) is equivalent to the 
10 minimization of E in equation (2a). In a similar manner, the secondary condition defined in 
equation (1) is equivalent to equations (13) and (17). The solution for equations (1 1)-(17) that 
can be determined with standard methods of numerical mathematics thereby simultaneously 
supplies the solution to the minimization problem of equations (1) and (2a), provided the 
measurement data is sufficiently great to ensure the secondary condition tj > 0. The latter can, 
1 5 however, be achieved in a simple manner by the addition of a constant to all measurement 
data. 

It is pointed out that there is a so-called permitted base vector, i.e. a vector a , 
which fulfils all boundary conditions. With the definition: 

r = max{m, |i=l,..JV} (19) 
20 the vector a may be composed of the following components: 

? = (r,..r) r (20) 
£ + =(0,...0) r (21) 
S~ ={r-m x ,,..T-m N ) T (22) 

A = (0,...0) r . (23) 
25 The method explained above may, in addition, be supplemented, in the context 

of a regularization, by a criterion that demands a smooth progression of the solution. In this 
regard, in addition to equation (2a), the sum of the absolute values of the second derivatives 
of the bolus arrival times t t is to be minimized below. 

To this end, a further set of variables S^.,.S^_ M _ X and S^ y ...S^_ M _ x is defined. 

30 Here, M stands for the number of vessel sections without successors, i.e. the downstream 
ends of the vessel tree under consideration. Notwithstanding the general principles, an 
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indexing in which these M vessel sections have the very indices N-M, . . . N-l may hereby be 
required. For each index i=l , ... (N-M-l), ST - 87 is to represent the second derivative of the 
bolus arrival times in the case of vessel section i. These variables therefore have to be linked 
to the vector of the bolus arrival times t by the following boundary condition: 

5 S:-S7=t m -2t t +|^r£'; with J = {j\p(j) = i} (24) 

Here, |J| corresponds to the number of vessel sections that have the vessel 
section i as their predecessor. Apart from junction points in the vessel tree, the vessel section 
i has only one successor j, so |J| = 1 and the right-hand side of equation (24) collapses for the 
usual discrete approximation of the second derivative. If the length of the individual vessel 
10 sections within the vessel tree varies greatly, the introduction of geometric weightings may 
be advisable. 

Linear programming may now be used to minimize the following variable: 

Zfe ++ ^")= Ski (2b) 
1=1 1=1 

A small value of E m hereby means that the absolute values of the second 
15 derivatives of the bolus arrival times are small. Therefore, "smooth" solutions t are 
preferred. 

The regularization further requires the introduction of a regularizing parameter 
X, which weights the two target functions to be minimized, i.e. the "data term" E in 
accordance with equation (2a) and the "model term" E m in accordance with equation (2b). 
20 The variables used for linear programming may then be assumed in 

accordance with equations (26) to (29) of Fig. 3, wherein I=N-M-1 is written for 
simplification purposes. The additional IxN coupling matrix C implements the right-hand 
side of the condition of equation (24) and is defined in accordance with equation (30) of Fig. 
3. 

25 During regularization, it is generally difficult to select a "good" value for the 

regularization parameter X in advance. To overcome this problem, a two-stage method is 
proposed here. In the first stage, the optimization is executed without regularization. This 
provides a first estimation ? (0) for the bolus arrival times, with which the "best possible" data 
term E (0> is then calculated, i.e. 

30 £ (O) = £h-', (0) | • (31) 
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A value for E that is achievable more practically may now be assumed, e.g. 
E = sE (0) with s > 1 . Since everything is linear, the original minimization of E may be 
replaced by the condition E = sE (0) , which leads to equations (32) to (35) in Fig. 3 for the 
linear programming. 

In accordance with equation (36), it is further possible to find a permitted base 
vector. It may be necessary, in turn, to shift all values ti (,) if the obtained t N (l) is negative. 

Alternatively, the worst case model term E m (0) may also be calculated, based 
on the first vector F <0) and the condition E m = sE m (0) with s < 1 may be added to the 
minimization problem of the data term E. 

Finally, it should be pointed out that the method described also permits the 
integration of other physiological information, which may be described in the form of linear 
boundary conditions. For example, the mass preservation could be formulated in this manner. 



